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Abstract: The normalization constant of the lateral distribution function (LDF) of an extensive air shower 
is a monotonous (almost linear) increasing function of the energy of the primary. Therefore, the interpo- 
lated signal at some fixed distance from the core can be calibrated to estimate the energy of the shower. 



There is, somehow surprisingly, a reconstructed optimal distance, 



opt 



at which the effects on the inferred 



signal, S(r op t), of the uncertainties on true core location, LDF functional form and shower-to-shower 
fluctuations are minimized. We calculate the value of r opt as a function of surface detector separation, 
energy and zenith angle and we demonstrate the advantage of using the r opt 



value of each individual 

shower instead of a same fixed distance for every shower, specially in dealing with events with saturated 
stations. The effects on the determined spectrum are also shown. 



Introduction 

In order to determine the energy of cosmic rays 
with surface detectors arrays, first the lateral dis- 
tribution function (LDF) of the shower particles, 
i.e. the particles density or signal versus distance 
to shower core location, is fitted assuming a known 
functional form. Following Hillas [1] proposal, the 
signal at some fixed distance of the shower core 
S(r op t) for all the showers, independent of their 
energy or direction, is used to relate it with primary 
energy, usually using monte carlo simulations. The 
optimum distance r opt is mainly related to the ge- 
ometry of the array. 

We show that this method may not reconstruct 
properly the shape of the spectrum. We use an 
AGASA-like experiment [2] as case study and in- 
ject both a single power law and a simplified, yet 
realistic, structured spectrum above 10 17 ' 7 eV. A 
special analysis has been done for saturated events. 



Method 

A previous version of our algorithm was presented 
in [3]. We use a simplified numerical approach to 
the simulation of extensive air shower detection in 
a surface array. Our ideal detector is an infinite 
array of equally spaced stations distributed in ele- 
mentary triangular cells of 1 km side. 

The input spectrum is a perfect power law spec- 
trum with index -3.0 from 10 18 to 10 20 - 1 and with 
an isotropic zenith distribution from to 45 de- 
grees. The number of events is 230.000, approx- 
imately the same statistic as the spectrum reported 
by AGASA. In a second case, we used a struc- 
tured spectrum (with second knee, ankle and GZK- 
cut-off, and exposure-limited at low energy) from 
10 17 - 7 to 10 210 eV as input. 

Firstly, we select random core position inside an 
elementary cell, and the signal of each station is 
estimated using the LDF reported by Auger [4]: 
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where /3(6>) = 3.1 — O.7sec(0). The signal ex- 
pected at each station is then fluctuated by Pois- 
sonian noise and recorded if it is above a threshold 
of 3.0 VEM (Vertical Equivalent Muons, the signal 
deposited by one vertical muon in an Auger water 
Cerenkov tank). Stations with a signal Sj > 5(0.2 
km , 1 EeV , 0°) are considered saturated and are 
excluded. 

The real core position is now moved using a gaus- 
sian distribution centered at this point. The sigma 
of this distribution is set taking into account the 
uncertainty in core determination, which depends 
on the array geometry and primary energy [5], The 
new core position simulates the reconstructed core. 
In order to mimic the reconstruction procedure we 
fit the signals of the triggered stations with the LDF 
of the AGASA form: 



log 5(r„ 



ai - a 2 log(r/r M ) - 0.61og(l + r 2 ) 

(2) 

with r&f = 91.6 m, the moliere radius at AGASA 
altitude. Finally, the signal at 600 m (5 e (600)), is 
used to estimate the primary energy. 

For the same shower, another procedure is done to 
obtained its real optimum distance r opt . The re- 
constructed core is shifted 50 times using again a 
gaussian distribution centered at this point and with 
the same sigma as before. For each new core posi- 
tion, an LDF is fitted and the point r opt is defined 
as the location of the minimum dispersion. With 
the signal at this point S(r opt ), the energy of the 
primary is again estimated. 



Conversion between Auger and AGASA 
LDFs 

A conversion of units between the LDF from 
Auger (that we have used to assign the signal in 
each station) and the LDF from AGASA (that we 
will use for energy determination) is needed. The 
AGASA LDF is: 

( - ^ - 1 ' 2 
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Figure 1 : r opt for different array spacings and pri- 
mary energies at 9 = 30 deg. The error bars repre- 
sent the CL at 68 and 95%. 



r 



-(r,(0)-1.2) 



(l + r 2 ) 



-0.6 



(3) 



where 77 = 3.84 - 2. 15(560(6*) - 1) and K is 
the shower size. The conversion factor (AGASA- 
LDF/Auger-LDF), depends on energy, zenith and 
core distance. A study of the energy and zenith de- 
pendence over the hole spectral range shows that it 
is negligible. Nevertheless, the dependence on core 
distance is sizable and a fit of the form l/(a+ bx c ) 
is used. 



' opt 



dependence with array spacing 



The dependence of r opt with the energy, zenith 
and detectors separation has been studied and it is 
presented in [6]. Here we show the dependence 
with array spacing for several primary energies and 
9 = 30 deg (see Fig. 1). The results are in agree- 
ment with the values used by AGASA (detectors 
separated 1 km and r opt = 600 m) and Auger (sep- 
aration of 1.5 km and r opt — 1000 m). Note, how- 
ever, the dependence of r opt with energy and its 
considerable dispersion. 

Energy determination 

Since traversed atmosphere is a function of zenith 
angle, AGASA experiment transformed the ob- 
served shower density 5^(600) at zenith angle 9 
into 5o(600), the corresponding value of a vertical 
shower. The attenuation curve is: 
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x exp 



S*(600) = S (600) / s (0) = S o (600) 
Ai A 2 



(4) 



where X = 920 g/cm 2 , Ai = 500 g/cm 2 and 
A2 = 594 g/cm 2 for showers with < 45 deg. 
The uncertainty in So (600) due to this transfor- 
mation is estimated to be ±5%. The conversion 
formula to relate So (600) with energy reported by 
AGASA is: 



E = 2.21 • 10 17 S (600) 1 U 'W 



(5) 



where different hadronic interaction models and 
simulation codes were considered. Using eq. 5 
we calculated the shower energy based on the ob- 
served signal at 600 m. 

In order to find the energy using the signal at r op t, 
we use the following parametrization of the shower 
size (obtained from eq. 3 and eq. 5): 

K(9, E EeV ) = 49.676 x f s (6) 

xh + ^°j X £1/1-03 (6) 

Once the S(r opt ) has been determined as explained 
before, the shower size is obtained from eq. 3 and 
the energy from eq. 6. 



Results and discussion 

First, using a flat spectrum from 10 17 7 to 10 210 
eV we calculated the distribution function of 
events contributing to each reconstructed energy 
both, for r(600) and r opt . The corresponding 68% 
and 95% confidence levels (at the low (L) and high 
energy sides (H)) for each distribution are shown 
in fig. 2.a-f. These distribution are very nearly 
Gaussian for r opt (fig.2.a-b) but skewed for r(600) 
(fig.2.c-d). This behaviour is somehow improved 
if events with saturated stations are eliminated 
(fig.2.e-f), although at the high cost of severely de- 
creasing high energy statistics (fig.2.h). Finally, we 
compare the median of the distributions with the 
corresponding reconstructed energy in order to as- 
sess the bias in each case (see, fig.2.g). In all cases 
the bias is negligible. To avoid border effects, last 
bins of the spectrum has been erased in the figures. 
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Figure 2: From a) to f): Fist row is for E(r op t), 
the second row is for E(600) with saturated events 
and the third row is for E(600) without saturated 
events; left column are the symmetry of the distri- 
butions comparing their high and low sides (Gaus- 
sian 0^ /ctl = 1) and right column are the tails of 
the distributions comparing their C.L. at 68% and 
95% (Gaussian 095/(768 = 2). (g) reconstructed 
energy bias, (h) fraction of saturated events. 



Figure 3 shows a reconstructed power law spec- 
trum from 10 18 to 10 20 1 . The slope of the spec- 
trum is better reconstructed using S(r opt ) than 
S(600). Again a considerable improvement is ob- 
tained by neglecting events with saturated stations, 
but at a high statistical cost. Furthermore, us- 
ing So (600) around 11% of the events are recon- 
structed outside of the input bounds (most of them 
are events in the lower energy bins). However, in 
the case of S(r op t) this is reduced to ~ 1% of the 
events. 

It is important to emphasize two things. First, in 
both cases, S#(600) and S(r opt ), the x 2 /ndof of 
the corresponding fits are very good (at the level 
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Figure 3: Spectrum input (black) with the same 
statistic as AGASA and the reconstructed one ob- 
tained using 5o(600) (S(r opt )) as energy estima- 
tor in blue (red). In green is the spectrum obtained 
using 5o(600) but rejecting saturated events. The 
vertical axis is multiplied by E 3 . Error bars are the 
CL at 68 and 95% with 1 12 spectrums. 
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Figure 4: Input structured spectrum (black) and 
the reconstructed spectra using 5*0(600) (blue) and 
S(r op t) (red) as energy estimator note the position 
of the ankle and the GZK transition. 
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of 10" 2 ). Second, in the case of 5 e (600), the 
energy reconstruction is very bad for events with 
one or more saturated stations and energy below 
10 18 ' 5 eV so they were rejected to improve re- 
construction. This problem does not happen with 
S(r op t), highlighting a major advantage of the lat- 
ter method. 

Figure 4 shows the effects of both energy recon- 
structions when applied to a realistically structured 
spectrum (black curve), that has a an ankle, a GZK 
modulation, observed by an array with tanh ac- 
ceptance that attains full efficiency < 10 18 eV. It 
can be seen that while the S(r opt ) method fairly 
reproduces the impinging spectrum, the Sg(600) 
method changes the position of the ankle and 
smoothes the GZK transition. 

Therefore, using an optimum distance, calculated 
for each individual shower to estimate primary en- 
ergy is a simple procedure to improve the reli- 
ability of the calculated spectrum. Additionally, 
the r opt strategy minimizes the dead-time intro- 
duced by saturated stations or the possible biases 
originated by the implementation of algorithms de- 
signed to recover saturated signals. 
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